Analisis de Datos Omicos en COVID
Alex Sanchez
Enero 2022
1 Introducción
1.1 Infraestructura informática para el análisis
Este estudio se llevará a cabo usando R/Bioconductor por lo que es preciso tener instaladas un conjunto de librerías. Esto puede hacerse siguiendo el procedimiento descrito a continuación.
El código que se presenta debería ejecutarse tan sólo una vez!
if (!require(BiocManager)) {
install.packages("BiocManager", dep = TRUE)
}
installifnot <- function(pckgName, BioC = TRUE) {
if (BioC) {
if (!require(pckgName, character.only = TRUE)) {
BiocManager::install(pckgName)
}
} else {
if (!require(pckgName, character.only = TRUE)) {
install.packages(pckgName, dep = TRUE)
}
}
}
installifnot("limma")
installifnot("edgeR")
installifnot("org.Hs.eg.db")
installifnot("clusterProfiler")
installifnot("dplyr", BioC = FALSE)
installifnot("gplots", BioC = FALSE)
installifnot("ggvenn", BioC = FALSE)
installifnot("pheatmap")
installifnot("prettydoc", BioC = FALSE)
installifnot("ggnewscale", BioC = FALSE)Puede resultar útil, aunque no es para nada imprescindible, disponer de, al menos dos, directorios específicos: - datos (o un nombre similar) donde guardar y de donde cargar los datos. - resultso (o un nombre similar) donde escribir los resultados.
2 Los datos para el análisis
2.1 Lectura de la matriz de contajes
Los datos de contajes se encuentran descargados desde GEO al archivo Rawcounts.csv desde donde se leeran.
Con la información sobre los grupos u otras covariables o variables auxiliares se creará un objeto (el habitual “targets”) que debe estar sincronizado con el anterior.
3 Preprocesado de los datos
3.1 Estandarización de los contajes
Además de filtrar, es bueno expresar los contajes en “CPMs” es decir “counts per million”, lo que no modificará los resultados del filtraje, pero estandarizará los valores, lo que es útil y necesario para los análisis posteriores.
COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0 0 0 0 0 0
ENSG00000227232 1 0 3 16 1 2
ENSG00000278267 3 1 2 0 3 0
ENSG00000243485 2 0 0 1 0 0
ENSG00000284332 0 0 0 0 0 0 COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973 0.0000000
ENSG00000243485 0.13289935 0.00000000 0.0000000 0.05301203 0.00000000 0.0000000
ENSG00000284332 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000Una vez los datos estan como CPMs, se procede a filtrarlos,
3.2 Filtraje de genes poco expresados
Los genes con recuentos muy bajos en todas las librerías proporcionan poca evidencia de expresión diferencial por lo que es habitual eliminar aquellos genes que, o bien son poco variables, o bien presentan poca o nula expresión en la mayoría de las muestras.
En este caso, siguiendo las indicaciones proporcionadas, se opta por conservar únicamente aquellos genes que presentan algún valor en, al menos, tres muestras de cada grupo.
[1] 60683 34[1] 33039 34Aunque no sea más que un ejemplo basta con ver los dos primeros genes para comprobar como el primero no cumple la condición, en el grupo “SANO”, mientras que los siguientes sí que la cumplen. Por lo tanto, al filtrar desaparece el primer gen de la matriz filtrada, pero no los dos siguientes.
COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
COV153 COV154 COV155 COV156 COV157 COV175
ENSG00000223972 0.000000 0.00000000 0.1163313 0.0000000 0.1173623 0.0000000
ENSG00000227232 1.388085 0.41957493 1.1633130 0.0000000 0.1173623 0.0000000
COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000223972 0.0000000 0.0000000 0.0575402 0.0000000 0.0000000 0.0000000
ENSG00000227232 0.9503052 1.9885554 0.1726206 0.1462825 0.6558384 0.3528731
HEA062 HEA063 HEA064 HEA065 HEA066 HEA067
ENSG00000223972 0.0000000 0.0000000 0.06183451 0.0000000 0.00000000 0.0000000
ENSG00000227232 0.0000000 0.0000000 0.06183451 0.1029819 0.05475139 0.0000000
HEA068 HEA069 HEA070 HEA071 HEA181 HEA182
ENSG00000223972 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
ENSG00000227232 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508 0.2800970
HEA183 HEA184 HEA185 HEA186
ENSG00000223972 0.00000000 0.000000 0.0000000 0.0000000
ENSG00000227232 0.05465108 0.000000 0.0000000 0.4151761
[ reached getOption("max.print") -- omitted 4 rows ] COV145 COV147 COV149 COV150 COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
COV152 COV153 COV154 COV155 COV156 COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
COV175 COV176 COV177 COV178 COV179 COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
HEA061 HEA062 HEA063 HEA064 HEA065 HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
HEA067 HEA068 HEA069 HEA070 HEA071
ENSG00000227232 0.00000000 0.06713275 0.06693886 0.04887949 0.1122282
ENSG00000278267 0.42685808 0.13426549 0.46857204 0.09775898 0.3366846
HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.20035083 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.20035083 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
[ reached getOption("max.print") -- omitted 4 rows ]3.3 Uso de clases específicas para manejar los datos
Cuando se trabaja con distintos objetos referidos a unos mismos datos, como la matriz de contajes y el objeto “targets”, es útil disponer de clases contenedoras que permitan trabajar con todos ellos a la vez, lo que no sólo facilita el trabajo sino que ayuda a evitar “desincronizaciones”.
Éste es el caso de la clase ExpressionSet habitualmente utilizada con microarrays o de la clase que la generaliza, llamada SummarizedExperiment.
Para datos de contaje es habitual usar una clase similar a ExpressionSet llamada DGEList” pensadas para manejar datos de contajes , definida en el paquete edgeR. Esta clase, más simple que las anteriores, utiliza listas para almacenar recuentos de “reads” e información asociada de tecnologías de secuenciación o expresión génica digital. Puede encontrarse información al respecto en la ayuda del paquete edgeR.
An object of class "DGEList"
$counts
COV145 COV147 COV149 COV150 COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
COV152 COV153 COV154 COV155 COV156 COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
COV175 COV176 COV177 COV178 COV179 COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
HEA061 HEA062 HEA063 HEA064 HEA065 HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
HEA067 HEA068 HEA069 HEA070 HEA071 HEA181
ENSG00000227232 0.0000000 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508
ENSG00000278267 0.4268581 0.13426549 0.46857204 0.09775898 0.3366846 0.2003508
HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$samples
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
29 more rows ...
$genes
genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...Uno de los aspectos interesantes de estas clases es la posibilidad de extraer partes de todos los objetos a la vez con el operador de “subsetting”.
[1] 33039 34 group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
COV157 COVID 999582.4 1 COV157 COVID red
COV175 COVID 998800.7 1 COV175 COVID red
COV176 COVID 999667.6 1 COV176 COVID red
COV177 COVID 999807.4 1 COV177 COVID red
COV178 COVID 999721.9 1 COV178 COVID red [1] "COV145" "COV147" "COV149" "COV150" "COV151" "COV157" "COV175" "COV176"
[9] "COV177" "COV178"Aunque podríamos haber creado el objeto a partir de todas las muestras, y haber realizado la extracción de genes y muestras posteriormente, hemos optado por no hacerlo para facilitar el seguimiento del proceso.
3.4 Normalización
Además de estandarizar los contajes, es importante eliminar otros los sesgos de composición entre librerías. Esto puede hacerse aplicando la normalización por el método TMM que genera un conjunto de factores de normalización, donde el producto de estos factores y los tamaños de librería definen el tamaño efectivo de la biblioteca.
La función calcNormFactors, de la librería edgeR, calcula los factores de normalización entre librerías.
Esto no modificará la matriz de contajes, pero actualizará los factores de normalización en el objeto DGEList (sus valores predeterminados son 1).
Es decir, aunque no se observen cambios en la matriz de contajes, cuando se utilizan estos factores de normalización en algún cálculo la importancia de las distintas columnas se tendrá en cuenta.
Resumiendo
Los análisis que se realicen a continuación se basaran en la matriz de contajes, filtrada, estandarizada y normalizada, sobre la que además se toman logaritmo base dos.
Esta será nuestra matriz de partida para los análisis siguientes,
4 Exploración de los datos
Una vez descartados los genes poco expresados y con los recuentos almacenados en un objeto DGEList, podemos`proceder a realizar algunos gráficos exploratorios para determinar si los datos aparentan buena calidad y/o si presentan algun problema.
4.1 Distribución de los contajes
Un diagrama de cajas con los datos, normalizados o no, muestra que la distribución de los contajes es muy asimétrica, lo que justifica la decisión de trabajar con los logaritmos de los datos.
La transformación logarítmica puede hacerse directamente pero es mejor usar la función cpm, como se ha hecho, que agrega una pequeña cantidad para evitar tomar logaritmos de cero.
4.2 Análisis de similaridad entre las muestras
4.2.1 Distancia entre muestras
La función dist permite calcular una matriz de distancias que contiene las comparaciones dos a dos entre todas las muestras. Por defecto se utiliza una distancia euclídea.
COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156
COV147 80.0
COV149 70.7 82.0
COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063
COV147
COV149
HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182
COV147
COV149
HEA183 HEA184 HEA185
COV147
COV149
[ reached getOption("max.print") -- omitted 31 rows ]Las matrices de distancias se pueden visualizar directamente con un heatmap.
Como puede verse _las muestras tienden a agruparse por el factor SANO/COVID, aunque una de las muestras. COV155 se separa del resto de las del grupo COVID.
4.2.2 Agrupamiento jerárquico
Un agrupamiento jerárquico proporciona una representación alternativa, también basada en la matriz de distancias.
Una de las muestras COVID parece más similar a las saludables que a las otras COVID.
4.2.3 Visualización en dimensión reducida
Como puede verse, el gráfico muestra la misma agrupación “natural” y el mismo comportamiento atípico de una muestra.
5 Análisis de expresión diferencial
El objetivo del análisis de expresión diferencial es seleccionar genes cuya expresión difiere entre grupos.
5.1 Selección de genes usando limma-Voom
La ventaja principal de esta aproximación es que permite trabajar con toda la flexibilidad de los modelos lineales para representar diseños experimentales, y, en muchos casos , aprovechar la experiencia previa del usuario en el manejo de limma.
5.1.1 Matriz de diseño y de contrastes
Utilizando la variable group podemos definir una matriz de diseño y, sobre ésta, los contrastes que nos interesan.
COVID SANO
COV145 1 0
COV147 1 0
COV149 1 0
COV150 1 0
COV151 1 0
COV152 1 0
COV153 1 0
COV154 1 0
COV155 1 0
COV156 1 0
COV157 1 0
COV175 1 0
COV176 1 0
COV177 1 0
COV178 1 0
COV179 1 0
COV180 1 0
HEA061 0 1
HEA062 0 1
HEA063 0 1
HEA064 0 1
HEA065 0 1
HEA066 0 1
HEA067 0 1
HEA068 0 1
HEA069 0 1
HEA070 0 1
HEA071 0 1
HEA181 0 1
HEA182 0 1
HEA183 0 1
HEA184 0 1
HEA185 0 1
HEA186 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"Dado que estamos interesados en las diferencias entre los grupos, necesitamos especificar qué comparaciones queremos llevar a cabo. Las comparaciones de interés se puede especificar utilizando la función makeContrasts. La matriz de contraste indica qué columnas de la matriz design vamos a comparar. En este caso tan sólo se llevará a cabo una comparación.
Contrasts
Levels CONTROLvsCOVID
COVID 1
SANO -15.1.2 Transformación de los contajes
An object of class "EList"
$genes
genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...
$targets
group lib.size norm.factors samples group.1 cols
COV145 COVID 1001653.3 1.0018278 COV145 COVID red
COV147 COVID 1100529.0 1.1010464 COV147 COVID red
COV149 COVID 955030.2 0.9553476 COV149 COVID red
COV150 COVID 947796.6 0.9479175 COV150 COVID red
COV151 COVID 899521.0 0.8996072 COV151 COVID red
29 more rows ...
$E
COV145 COV147 COV149 COV150 COV151
ENSG00000227232 -0.8223650 -1.1381984 -0.4037625 0.50837551 -0.6500950
ENSG00000278267 -0.5183002 -0.9448596 -0.5597126 -0.92265092 -0.3219038
COV152 COV153 COV154 COV155 COV156
ENSG00000227232 -0.6246522 0.9358596 0.1073289 1.0159241 -0.8282289
ENSG00000278267 -1.0045156 -0.9810644 -0.5828556 -0.2862517 -0.4994115
COV157 COV175 COV176 COV177 COV178
ENSG00000227232 -0.74055345 -0.9859431 0.5497283 1.3062240 -0.6882933
ENSG00000278267 -0.74055345 -0.4526424 -0.4619603 -0.1060517 -1.1161581
COV179 COV180 HEA061 HEA062 HEA063
ENSG00000227232 -0.654532686 0.1342634 -0.20656331 -1.01143480 -1.1012204
ENSG00000278267 -0.499855856 -1.0746764 -0.42066201 -0.35512422 -0.6309653
HEA064 HEA065 HEA066 HEA067 HEA068
ENSG00000227232 -0.9265423 -0.7259291 -0.78574256 -1.1218171 -0.9590727
ENSG00000278267 -0.1024608 -0.6077105 0.05386471 -0.2313968 -0.7976723
HEA069 HEA070 HEA071 HEA181 HEA182
ENSG00000227232 -1.0099326 -0.8215555 -0.6475041 -0.4772445 -0.4479095
ENSG00000278267 -0.2372664 -0.6984811 -0.1968898 -0.4772445 -0.7333818
HEA183 HEA184 HEA185 HEA186
ENSG00000227232 -0.8517698 -1.0495769 -0.9932084 -0.1628122
ENSG00000278267 -0.2738825 0.2117625 -0.3086203 -0.5957197
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$weights
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 6.976324 6.423290 7.345017 7.414329 7.913732 6.966953 7.070951
[2,] 13.344117 9.792171 16.233083 16.785098 21.291716 13.276959 14.091120
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 8.568175 9.129126 8.046130 6.792973 7.049157 7.046103 6.946915
[2,] 27.856177 34.175509 22.694541 12.084889 13.891449 13.863695 13.134382
[,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 6.497274 6.878669 6.667025 21.376991 19.139568 14.546448 14.827401
[2,] 10.280091 12.659254 11.283734 8.528972 8.192672 7.560903 7.603134
[,22] [,23] [,24] [,25] [,26] [,27] [,28]
[1,] 20.095913 24.685105 13.715056 13.110198 11.661751 22.943084 24.332783
[2,] 8.340050 8.955089 7.428218 7.308304 7.031769 8.740787 8.912791
[,29] [,30] [,31] [,32] [,33] [,34]
[1,] 22.360747 15.055081 19.758022 16.994284 20.284112 17.780471
[2,] 8.666115 7.636840 8.288627 7.906856 8.368394 8.008534
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$design
COVID SANO
COV145 1 0
COV147 1 0
COV149 1 0
COV150 1 0
COV151 1 0
29 more rows ...5.1.3 Selección de genes diferencialmente expresados
Como en el caso de los microarrays el objeto voomObj y las matrices de diseño y contrastes se utilizaran para ajustar un modelo y, a continuación realizar las comparaciones especificadas sobre el modelo ajustado. El proceso finaliza con la regularización del estimador del error usando la función eBayes.
5.1.4 Top tables
Los resultados de un análisis de expresión diferencial se pueden extraer con la función topTable. Esta función genera una tabla de resultados cuyas columnas contienen información acerca de los genes y la diferencia entre los grupos comparados. Concretamente:
genes logFC AveExpr t P.Value
ENSG00000167900 ENSG00000167900 2.832423 4.111065 12.36669 8.088344e-15
ENSG00000090889 ENSG00000090889 2.476447 1.292066 11.83766 3.015984e-14
ENSG00000111206 ENSG00000111206 2.423923 2.018435 11.80083 3.309678e-14
ENSG00000166851 ENSG00000166851 3.223189 2.615540 11.65570 4.781203e-14
ENSG00000135476 ENSG00000135476 2.236975 1.705254 11.53985 6.425179e-14
ENSG00000123485 ENSG00000123485 2.553959 1.914888 11.46276 7.828950e-14
adj.P.Val B
ENSG00000167900 2.672308e-10 23.54910
ENSG00000090889 3.644948e-10 22.09352
ENSG00000111206 3.644948e-10 21.99541
ENSG00000166851 3.949154e-10 21.70990
ENSG00000135476 4.013294e-10 21.35292
ENSG00000123485 4.013294e-10 21.189535.1.5 Visualización de los resultados
Para visualizar los resultados podemos usar un volcanoPlot:
Con el fin de observar si existen perfiles de expresión diferenciados podemo realizar un mapa de colores con los genes más diferencialmente expresados.
Es decir, fijamos un criterio de selección de genes y retenemos aquellos componentes de la tabla de resultados que lo cumplen. Por ejemplo: Genes con un p-balor ajustado inferior a 0.001 y un `fold-change’ superior a 2.
[1] 199Con la matriz de expresión de los genes que verifican dicha condición se puede construir un heatmap.
Los dos grupos estan diferenciados, sobre todo en un subconjunto de genes en donde las expresiones toman signos (colores) distintos. En otro de los grupos parece que, en el grupo de individuos sanos los genes diferencialmente expresados se encuentran sobre-expresados en el grupo COVID, y apenas expresados en el grupo sanos.
5.2 Análisis de expresión diferencial usando el paquete edgeR
El análisis con edgeR es similar al anterior (se originan en el mismo equipo de investigación) pero la modelización es distinta.
El análisis utiliza un GLM pero, en una forma que recuerda a lo que se hace con limma, realiza un paso adicional,en el que se calcula una estimación mejorada de la dispersión (variabilidad) de las muestras que integra las estimaciones individuales y la global mediante estimación Bayes empírica.
Con este objeto, que añade a los contajes normalizados los estimadores mejorados de dispersión, se ajusta un modelo lineal generalizado con distribución binomial para los errores.
Una vez ajustado el modelo se procede a construir el contraste y realizar el test.
De hecho podemos usar la matriz de contrastes que construímos para limma voom, en la misma forma que hemos reutilizado la de diseño.
An object of class "DGELRT"
$coefficients
COVID SANO
ENSG00000227232 -14.28149 -15.28977
ENSG00000278267 -15.16817 -14.68027
ENSG00000233750 -15.61279 -15.65247
ENSG00000268903 -13.79667 -14.19322
ENSG00000269981 -14.09073 -14.31982
ENSG00000241860 -14.96062 -14.97429
$fitted.values
COV145 COV147 COV149 COV150 COV151
ENSG00000227232 0.50362351 0.55333742 0.48018176 0.47654477 0.45227218
ENSG00000278267 0.13404573 0.14727771 0.12780641 0.12683838 0.12037793
COV152 COV153 COV154 COV155 COV156
ENSG00000227232 0.50436793 0.49623562 0.42920617 0.41355987 0.4463544
ENSG00000278267 0.13424387 0.13207935 0.11423862 0.11007416 0.1188028
COV157 COV175 COV176 COV177 COV178
ENSG00000227232 0.51862932 0.49791657 0.49815307 0.50596775 0.5449482
ENSG00000278267 0.13803971 0.13252676 0.13258970 0.13466968 0.1450448
COV179 COV180 HEA061 HEA062 HEA063
ENSG00000227232 0.51149863 0.52950246 0.10257076 0.1050509 0.11179640
ENSG00000278267 0.13614179 0.14093373 0.29173133 0.2987852 0.31797088
HEA064 HEA065 HEA066 HEA067 HEA068
ENSG00000227232 0.11129685 0.10394130 0.09967537 0.11340392 0.11490841
ENSG00000278267 0.31655006 0.29562943 0.28349630 0.32254297 0.32682204
HEA069 HEA070 HEA071 HEA181 HEA182
ENSG00000227232 0.11899089 0.10109909 0.09995140 0.10161035 0.11090222
ENSG00000278267 0.33843341 0.28754563 0.28428137 0.28899975 0.31542764
HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.10432431 0.10786525 0.10373204 0.10677595
ENSG00000278267 0.29671878 0.30678991 0.29503426 0.30369172
[ reached getOption("max.print") -- omitted 4 rows ]
$deviance
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981
13.203208 3.466706 1.937466 6.511782 10.312153
ENSG00000241860
4.314431
$method
[1] "oneway"
$unshrunk.coefficients
COVID SANO
ENSG00000227232 -14.50309 -16.07675
ENSG00000278267 -15.82674 -15.03147
ENSG00000233750 -17.00917 -17.18018
ENSG00000268903 -13.92723 -14.39407
ENSG00000269981 -14.27015 -14.55116
ENSG00000241860 -15.45818 -15.48078
$df.residual
[1] 32 32 32 32 32 32
$design
COVID SANO
COV145 1 0
COV147 1 0
COV149 1 0
COV150 1 0
COV151 1 0
29 more rows ...
$offset
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 13.81716 13.9113 13.7695 13.7619 13.70962 13.81864 13.80238 13.65727
[,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,] 13.62014 13.69645 13.84652 13.80577 13.80624 13.82181 13.89602 13.83268
[,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
[1,] 13.86727 13.79954 13.82344 13.88567 13.88119 13.81282 13.77091 13.89995
[,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
[1,] 13.91313 13.94804 13.78509 13.77367 13.79014 13.87764 13.8165 13.84987
[,33] [,34]
[1,] 13.8108 13.83972
attr(,"class")
[1] "CompressedMatrix"
attr(,"Dims")
[1] 6 34
attr(,"repeat.row")
[1] TRUE
attr(,"repeat.col")
[1] FALSE
$dispersion
[1] 9.765625e-05 9.765625e-05 9.765625e-05 9.765625e-05 9.765625e-05
[6] 9.765625e-05
$prior.count
[1] 0.125
$AveLogCPM
[1] 1.198866 1.145721 1.023866 1.444349 1.352131 1.128949
$df.residual.zeros
[1] 32 32 32 32 32 32
$df.prior
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981
6.412772 9.221763 9.221763 9.221763 9.221763
ENSG00000241860
9.221763
$var.post
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981
0.3701341 0.1170509 0.0646997 0.2054864 0.2930194
ENSG00000241860
0.1356745
$var.prior
[1] 0.15822624 0.14729701 0.07911392 0.21240299 0.19157096 0.13861893
$samples
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1.0018278 COV145 COVID red
COV147 COVID 999530.1 1.1010464 COV147 COVID red
COV149 COVID 999667.8 0.9553476 COV149 COVID red
COV150 COVID 999872.5 0.9479175 COV150 COVID red
COV151 COVID 999904.2 0.8996072 COV151 COVID red
29 more rows ...
$genes
genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
ENSG00000241860 ENSG00000241860
$table
logFC logCPM F PValue
ENSG00000227232 1.45463439 1.198866 13.174669854 0.02722442
ENSG00000278267 -0.70389442 1.145721 9.137260461 0.30104772
ENSG00000233750 0.05724585 1.023866 0.145023951 0.92283279
ENSG00000268903 0.57210237 1.444349 6.405756530 0.25124362
ENSG00000269981 0.33050524 1.352131 1.266644864 0.54237109
ENSG00000241860 0.01973312 1.128949 0.006129263 0.97699496
$comparison
[1] "1*COVID -1*SANO"
$df.test
[1] 1 1 1 1 1 1
$df.total
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981
38.41277 41.22176 41.22176 41.22176 41.22176
ENSG00000241860
41.22176 Los resultados se almacenan en un objeto similar a la ’ topTable’ de limma.
Coefficient: 1*COVID -1*SANO
genes logFC logCPM F PValue
ENSG00000123485 ENSG00000123485 2.969009 2.874508 143.3978 6.410064e-15
ENSG00000167900 ENSG00000167900 3.052165 4.946915 165.9368 1.024164e-14
ENSG00000187837 ENSG00000187837 1.719557 4.813853 135.1654 1.889383e-14
ENSG00000111206 ENSG00000111206 2.834520 2.919472 132.0747 2.767000e-14
ENSG00000178999 ENSG00000178999 2.987952 3.225967 135.0640 4.056862e-14
ENSG00000166851 ENSG00000166851 3.638997 3.745942 144.7154 5.207210e-14
FDR
ENSG00000123485 1.691868e-10
ENSG00000167900 1.691868e-10
ENSG00000187837 2.080777e-10
ENSG00000111206 2.285473e-10
ENSG00000178999 2.680693e-10
ENSG00000166851 2.867350e-10Podemos seleccionar los genes más diferencialmente expresados de la misma forma que hicimos con limma-voom
[1] 376Obsérvese que la lista de genes seleccionados es muy similar en ambos casos.
6 Anotación de resultados y análisis de significación biológica
Para el análisis de significación se utilizan dos listas de transcritos:
- La lista de transcritos diferencialmenete expresados
- La lista de todos los tránscritos o “Universo”
[1] 383[1] 330396.1 Anotación de los identificadores
Esto es posible, y de hecho sencillo de llevar a cabo, usando el paquete annotate.
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GENETYPE" "GO" "GOALL" "IPI" "MAP"
[16] "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM"
[21] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
[26] "UNIPROT" ENSEMBL SYMBOL ENTREZID GENENAME
1 ENSG00000167900 TK1 7083 thymidine kinase 1
2 ENSG00000090889 KIF4A 24137 kinesin family member 4A
3 ENSG00000111206 FOXM1 2305 forkhead box M1
4 ENSG00000166851 PLK1 5347 polo like kinase 1
5 ENSG00000135476 ESPL1 9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485 HJURP 55355 Holliday junction recognition protein[1] 385 4Como puede verse, el número de anotaciones es el mismo que el de identificadores ENSEMBL, lo que podría llevar a pensar que, es posible que, antes de subir los datos a GEO se hayan agrupado los contajes por genes.
Para la anotación del universo se procederá igual.
ENSEMBL SYMBOL ENTREZID GENENAME
1 ENSG00000167900 TK1 7083 thymidine kinase 1
2 ENSG00000090889 KIF4A 24137 kinesin family member 4A
3 ENSG00000111206 FOXM1 2305 forkhead box M1
4 ENSG00000166851 PLK1 5347 polo like kinase 1
5 ENSG00000135476 ESPL1 9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485 HJURP 55355 Holliday junction recognition protein[1] 33202 4En este caso se observa como hay más anotaciones que transcritos, lo que sugiere que múltiples transcritos han sido mapeados en el mismo gen.
6.2 Análisis de enriquecimiento
El paquete clusterProfiler admite identificadores de tipo ENSEMBL y permite gran variedad de análisis complementarios al enriquecimiento por lo que, es una de las mejores opciones para el análisis de significación biológica.
ID
GO:0006958 GO:0006958
GO:0002455 GO:0002455
GO:0006956 GO:0006956
GO:0016064 GO:0016064
GO:0019724 GO:0019724
Description
GO:0006958 complement activation, classical pathway
GO:0002455 humoral immune response mediated by circulating immunoglobulin
GO:0006956 complement activation
GO:0016064 immunoglobulin mediated immune response
GO:0019724 B cell mediated immunity
GeneRatio BgRatio pvalue p.adjust qvalue Count
GO:0006958 80/342 120/16130 4.046401e-107 1.215134e-103 1.117659e-103 80
GO:0002455 81/342 135/16130 9.256747e-103 1.389901e-99 1.278405e-99 81
GO:0006956 81/342 148/16130 4.506649e-98 4.511156e-95 4.149280e-95 81
GO:0016064 81/342 200/16130 3.900362e-84 2.928197e-81 2.693303e-81 81
GO:0019724 81/342 203/16130 1.728105e-83 1.037900e-80 9.546418e-81 81Con los resultados del análisis de enriquecimiento se pueden llevar a cabo distintas visualizaciones cuya interpretación exacta puede verse en el manual de clusterProfiler
